function [R_maneuver,V_maneuver,t_maneuver,delV_total] = Run_StationKeeping(S_chase_eci,S_target_eci,S_target_ecef,...
                                        t_crossing,MAX_STATIONKEEPING_DURATION_HOUR,UTC_mission_start);
mu = 398600.4418;
%% Find the states of the chase vehicle in HCW frame of the S_target_eci;
% Define new variables

Rt = S_target_eci(1:3);
Vt = S_target_eci(4:6);
Rc = S_chase_eci(1:3);
Vc = S_chase_eci(4:6);

% Define HCW frame

Ht = cross(Rt,Vt);
%[km^2/s]Target vehicle specific angular momentum.

Wt = Ht / dot(Rt,Rt);
%[rad/s]Target vehicle angular velocity.

Rhat = Rt / norm(Rt);
%[]Current radial direction.

Nhat = cross(Rt,Vt) / norm(cross(Rt,Vt));
%[]Current normal direction.

That = cross(Nhat,Rhat);
%[]Current tangential direction.

Hcw2Eci = [Rhat, That, Nhat];
%[]Matrix that transforms vectors from HCW coordinates to ECI coordinates.

Rct_rel = transpose(Hcw2Eci) * (Rc - Rt);
%[km]Current relative position in HCW coordinates.

Wt = transpose(Hcw2Eci) * Wt;
%[rad/s]Transforms the target vehicle's rotational velocity from ECI to HCW coordinates.

Vct_rel = transpose(Hcw2Eci) * (Vc - Vt - cross(Wt,Rct_rel));
%[km/s]Current relative velocity in HCW coordinates.

%% Determine best maneuver between given time interval

MAX_STATIONKEEPING_DURATION_SEC = MAX_STATIONKEEPING_DURATION_HOUR*3600;

n = sqrt(mu/norm(Rt)^3);

deltat = 3600 : 100 : MAX_STATIONKEEPING_DURATION_SEC;
totaldeltaV = zeros(1,length(deltat));
deltaV1 = zeros(3,length(deltat));
deltaV2 = zeros(3,length(deltat));
for k=1:length(deltat)
    
    F = [4-3*cos(n*deltat(k)) , 0 , 0;...
        6*(sin(n*deltat(k)) - n*deltat(k)) , 1 , 0;...
        0 , 0 , cos(n*deltat(k))];
    %[] F matrix
    
    G = [(1/n)*sin(n*deltat(k)) , (2/n)*(1-cos(n*deltat(k))) , 0;...
        (2/n)*(cos(n*deltat(k)) - 1) , (1/n)*(4*sin(n*deltat(k)) - 3*n*deltat(k)) , 0;...
        0 , 0 , (1/n)*sin(n*deltat(k))];
    %[] G matrix
    
    Vct_rel_firstburn = G \ (-F * Rct_rel);
    %[] Velocity after the first burn
    
    deltaV1(:,k) = Vct_rel_firstburn - Vct_rel;
    %[] change in velocity 1
    
    Fdot = [3*n*sin(n*deltat(k)) , 0 , 0;...
        6*(n*cos(n*deltat(k))-n) , 0 , 0;...
        0 , 0 , -n*sin(n*deltat(k))];
    %[] F dot matrix. Derivative respect to deltat(k)
    
    Gdot = [cos(n*deltat(k)) , 2*sin(n*deltat(k)) , 0;...
        -2*sin(n*deltat(k)) , 4*cos(n*deltat(k))-3 , 0 ;...
        0 , 0 , cos(n*deltat(k))];
    %[] G dot matrix. Derivative respect to deltat(k)
    
    Vct_rel_secondburn = Fdot*Rct_rel + Gdot*(Vct_rel_firstburn);
    %[] Second burn at the meeting
    
    deltaV2(:,k) = -Vct_rel_secondburn;
    %[] delta V required for the second burn
    
    totaldeltaV(k) = abs(norm(deltaV1(:,k))) + abs(norm(deltaV2(:,k)));
end

delV = min(totaldeltaV);

min_index = find(totaldeltaV==delV);

Vct_rel_plus = Vct_rel+deltaV1(:,min_index);

deltat_min = deltat(min_index);

Vc_plus = Hcw2Eci*Vct_rel_plus+Vt+cross(Wt,Rct_rel);

Sc_plus = [Rc;Vc_plus];

Sc_plus_meters = Sc_plus * 1000;

S_maneuver = Ephemeris_new(Sc_plus_meters,t_crossing,t_crossing+deltat_min);
S_maneuver = S_maneuver';
t = S_maneuver(1,:);
R = S_maneuver(2:4,:)/1000;
V = S_maneuver(5:7,:)/1000;
% this includes states of the spacecraft after maneuver
%% Find target's state after maneuver

UTC_tf = [UTC_mission_start(1:5),UTC_mission_start(6)+t(end)];
year = UTC_tf(1);
mon = UTC_tf(2);
day = UTC_tf(3);
hour = UTC_tf(4);
minute = UTC_tf(5);
sec = UTC_tf(6);

Mjd_UTC_tf = Mjday(year, mon, day, hour, minute, sec);

S_target_eci_tf = ECEF2ECI(Mjd_UTC_tf,S_target_ecef')/1000;

Rt_tf = S_target_eci_tf(1:3);
Vt_tf = S_target_eci_tf(4:6);
Rc_tf = R(:,end);
Vc_tf = V(:,end);
% Define HCW frame

Ht = cross(Rt_tf,Vt_tf);
%[km^2/s]Target vehicle specific angular momentum.

Wt = Ht / dot(Rt_tf,Rt_tf);
%[rad/s]Target vehicle angular velocity.

Rhat = Rt_tf / norm(Rt_tf);
%[]Current radial direction.

Nhat = cross(Rt_tf,Vt_tf) / norm(cross(Rt_tf,Vt_tf));
%[]Current normal direction.

That = cross(Nhat,Rhat);
%[]Current tangential direction.

Hcw2Eci_tf = [Rhat, That, Nhat];
%[]Matrix that transforms vectors from HCW coordinates to ECI coordinates.

Rct_rel_tf = transpose(Hcw2Eci_tf) * (Rc_tf - Rt_tf);
%[km]Current relative position in HCW coordinates.

Wt = transpose(Hcw2Eci_tf) * Wt;
%[rad/s]Transforms the target vehicle's rotational velocity from ECI to HCW coordinates.

Vct_rel = transpose(Hcw2Eci_tf) * (Vc_tf - Vt_tf - cross(Wt,Rct_rel_tf));
%[km/s]Current relative velocity in HCW coordinates.

Vct_rel_plus_tf = Vct_rel+deltaV2(:,min_index);

Vc_plus_tf = Hcw2Eci_tf*Vct_rel_plus_tf+Vt_tf+cross(Wt,Rct_rel_tf);

S_plus_tf = [Rc_tf;Vc_plus_tf];

delV1_eci = norm(Vc-Vc_plus);

delV2_eci = norm(Vc_tf-Vc_plus_tf);

%% save output
numb_data = length(t);
R_maneuver = [R(:,1:numb_data-1),S_plus_tf(1:3)];
V_maneuver = [V(:,1:numb_data-1),S_plus_tf(4:6)];
t_maneuver = t;
delV_total = [delV1_eci,delV2_eci];























end